Quantile Deviation Efficiency Visualization

Libraries, Functions, and Data

source('../ttm_to_map/functions.R')

library(tidyverse)
library(leaflet)
library(sf)
library(glue)
library(reshape2)
library(mapview)


original_scores_long <- read.csv('../../data/html_mapping_data/scores_long.csv')

Compute Efficiency (Deviation)

Efficiency: the discrepancy between population and accessibility percentiles.

The larger the discrepancy, the less efficient the transit system is in a particular block. This assumes that maximal efficiency maximizes transit resources where there are higher populations, and vice versa.

Disadvantage: locations that people commonly access using transit without a large population are penalized when they shouldn’t be. For example schools and high density amenity regions. To improve the model we need to consider not only population but also visitor traffic in each block.

# filter to keep only unweighted and nearest 2 scores
scores <- original_scores_long %>% subset(weight == 'no' & nearest_n == '2')

# function to convert vector to percentile
vec_to_percentile <- function(vec) {
  percentile <- ecdf(vec)
  return(round(percentile(vec), 4)*100)
}

# population percentiles
scores$pop_percentile <- vec_to_percentile(scores$pop)

# score percentiles for each amenity
scores <- scores %>% group_by(type) %>% summarise(fromId, type, score, score_percentile = vec_to_percentile(score), pop, pop_percentile) %>% arrange(fromId)
scores <- scores[c(2,1,3,4,5,6)]

# discrepancy is the efficiency
# linear penalization weighs all differences equally
scores$abs_deviation <- -1*round(abs(scores$score_percentile - scores$pop_percentile), 2)
# quadratic penalization weighs larger discrepencies more
scores$sqr_deviation <- -1*round((scores$score_percentile - scores$pop_percentile)^2, 2)

# keep only relevant columns
scores <- scores[-3]

# view deviation distribution
x <- scores %>%
  ggplot(aes(x = abs_deviation, color = type)) +
  geom_density() +
  egg::theme_article() +
  theme(aspect.ratio = 0.4) +
  ggtitle('Absolute Deviation Density')

y <- scores %>%
  ggplot(aes(x = sqr_deviation, color = type)) +
  geom_density() +
  egg::theme_article() +
  theme(aspect.ratio = 0.4) +
  ggtitle('Square Deviation Density')

gridExtra::grid.arrange(x, y)

scores_long <- melt(scores, id.vars=c(-6,-7), value.name = 'deviation', variable.name = 'deviation_type')

sample_n(scores_long, 5)

Bind Deviation to Shape Data

# import shapees
canada_shape <- st_read("../../data/census2016_DBS_shp/DB_Van_CMA/DB_Van_CMA.shp", stringsAsFactors = FALSE)

## Reading layer `DB_Van_CMA' from data source `C:\Users\Luka\Documents\GitHub\Capstone-Project\data\census2016_DBS_shp\DB_Van_CMA\DB_Van_CMA.shp' using driver `ESRI Shapefile'
## Simple feature collection with 15197 features and 27 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4001643 ymin: 1957237 xmax: 4068119 ymax: 2032663
## Projected CRS: PCS_Lambert_Conformal_Conic

# select a greater metropolitan area
metropolitan_area <- "Vancouver"

# filter columns and rows
vancouver_shape <- data.frame(canada_shape[which(canada_shape$CMANAME == metropolitan_area), c(1, 28)])

# id to factor
vancouver_shape$DBUID <- as.double(vancouver_shape$DBUID)

# join factor and geometry data 
scores_viz_frame <- left_join(vancouver_shape, scores_long, by = c('DBUID' = 'fromId'))

# convert back to sf object
scores_viz_frame_sf <- st_as_sf(scores_viz_frame)

# convert to st object
scores_viz_frame_st <- st_transform(scores_viz_frame_sf, crs = 4326)

Visualization Function

# takes an st data object, amenity type, deviation type and displays an efficiency map
# deviation = 1 is linear penalization
# deviation = 2 is square penalization
map_maker_efficiency <- function(data, amenity, deviation_type = 2, output_dir, view_map = TRUE) {
  
  amn_name <- amenity %>%
              str_to_title() %>%
              str_replace_all('Or', 'or') %>%
              str_replace('And', 'and') %>%
              str_replace('/Performance', '')

  deviation_type <- c('abs_deviation', 'sqr_deviation')[deviation_type]
  
  file_name <- glue('{amn_name} Transit Efficiency - Deviation ({substr(deviation_type, 1, 3)})')
  print(paste('Current Map:', file_name))

  # subset info
  polyg_subset <- data[data$type == amenity & data$deviation_type == deviation_type, ]
  
  # score vector
  deviation_vec <- polyg_subset$deviation
  
  # popup 
  dev_percentile <- ecdf(deviation_vec)
  p_popup <- paste0(
                    "Efficiency Percentile: <strong>", round(dev_percentile(deviation_vec), 2)*100, '%',"</strong>", "<br>",
                    "<i>The higher the efficiency percentile, the better.</i>", "<br><br>",
                    
                    "Deviation (", deviation_type, "):  <strong>", -1*polyg_subset$deviation, " </strong><br><br>",
                    
                    "Population Percentile:  <strong>", polyg_subset$pop_percentile, " </strong><br>",
                    "Access Score Percentile:  <strong>", polyg_subset$score_percentile, " </strong><br>",
                    "Block ID: ", polyg_subset$DBUID, "<br><br>",
                    
                    "Notes: <br>",
                    "<i>abs_deviation = absolute(population% - score%)</i>", "<br>",
                    "<i>sqr_deviation = (population% - score%)^2</i>")
                    
  # colour palette 
  Rd2Gn <- c("#e30606", "#fd8d3c", "#ffe669", "#cdff5e", "#64ed56")
  pal_fun <- colorQuantile(palette = Rd2Gn, NULL, n = 5)
  #Rd2Wh <- c("#ff1919", "#ff6666", "#ffb2b2", "#ffe5e5", "#ffffff")
  #pal_fun <- colorQuantile(palette = Rd2Wh, NULL, n = 5)
  
  map <- leaflet(data = polyg_subset) %>%
      addPolygons(
        stroke = FALSE, 
        fillColor = ~pal_fun(deviation_vec), 
        fillOpacity = 0.6, # smoothFactor = 0.0 --> smooth factor does nothing 
        popup = p_popup) %>%
      addTiles() %>%
      setView(lng = -122.8, lat = 49.2, zoom = 11) %>%
      addLegend("bottomleft",  
                pal = pal_fun,    
                values = ~deviation_vec,
                title = glue('{amn_name} Transit Efficiency'))
    
  if (view_map == TRUE) {
    return(map)
  } else {
    mapshot(map, url = glue("{getwd()}/{output_dir}/{file_name}.html"))
  }
  
}

map_maker_efficiency(scores_viz_frame_st,
                     amenity = "museum",
                     deviation = 1,
                     view_map = TRUE)

## [1] "Current Map: Museum Transit Efficiency - Deviation (abs)"

map_maker_efficiency(scores_viz_frame_st,
                     amenity = "museum",
                     deviation = 2,
                     view_map = TRUE)

## [1] "Current Map: Museum Transit Efficiency - Deviation (sqr)"

Export

types <- unique(scores_long$type)

for(deviation in 1:2) {
  for (type in 1:4) {
    map_maker_efficiency(scores_viz_frame_st, amenity = types[type], deviation = deviation,
                         view_map = FALSE, output_dir = '../../data/html_maps/quantile_deviation_efficiency_maps')
  }
}